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ABSTRACT 

Context. Recent simulations have explored different ways to form accretion disks around low-mass stars. However, it has been difficult 
to differentiate between the proposed mechanisms because of a lack of observable predictions from these numerical studies. 

Aims. We aim to present observables that can differentiate a rotationally supported disk from an infalling rotating envelope toward 
deeply embedded young stellar objects (Menv > ^disk) and infer their masses and sizes. 

Methods. Two 3D magnetohydrodynamics (MHD) formation simulations of Li and collaborators are studied, with a rotationally sup¬ 
ported disk (RSD) forming in one but not the other (where a pseudo-disk is formed instead), together with the 2D semi-analytical 
model. The dust temperature structure is determined through continuum radiative transfer RADMC3D modelling. A simple tempera¬ 
ture dependent CO abundance structure is adopted and synthetic spectrally resolved submm rotational molecular lines up to /u = 10 
are compared with existing data to provide predictions for future ALMA observations. 

Results. 3D MHD simulations and 2D semi-analytical model predict similar compact components in continuum if observed at the 
spatial resolutions of 0.5-1" (70-140 AU) typical of the observations to date. A spatial resolution of ~14 AU and high dynamic range 
(> 1000) are required in order to differentiate between RSD and pseudo-disk formation scenarios in the continuum. The moment 
one maps of the molecular lines show a blue- to red-shifted velocity gradient along the major axis of the flattened structure in the 
case of RSD formation, as expected, whereas it is along the minor axis in the case of a pseudo-disk. The peak-position velocity 
diagrams indicate that the pseudo-disk shows a flatter velocity profile with radius than an RSD. On larger-scales, the CO isotopolog 
line profiles within large (> 9") beams are similar and are narrower than the observed line widths of low-/ (2-1 and 3-2) lines, 
indicating significant turbulence in the large-scale envelopes. However a forming RSD can provide the observed line widths of high-/ 
(6-5, 9-8, and 10-9) lines. Thus, either RSDs are common or a higher level of turbulence (b ~ 0.8 km s“^) is required in the inner 
envelope compared with the outer part (0.4 km s"^). 

Conclusions. Multiple spatially and spectrally resolved molecular line observations can differentiate between the pseudo-disk and 
the RSD much better than continuum data. The continuum data give a better estimate on disk masses whereas the disk sizes can 
be estimated from the spatially resolved molecular lines observations. The general observable trends are similar between the 2D 
semi-analytical models and 3D MHD RSD simulations. 

Key words, stars: formation, radiative transfer, accretion disks; line: profiles; methods:numerical, magnetohydrodynamics (MHD) 


1. Introduction 

The formation of stars and their planetary systems is linked 
through the formation and evolution of accretion disks. In the 
standard star formation picture, the infalling material forms an 
accretion disk simply from angular momentum conservation 
(e.g., Lin & Pringle 1990; Bodenheimer 1995; Belloche 2013). 
However, magnetic field strengths observed toward molecular 
cores (see Crutcher 2012, for a recent review) are expected the¬ 
oretically to be sufficient in affecting the formation and evolu¬ 
tion of disks around low-mass stars (e.g., Galli et al. 2006; Joos 
et al. 2012; Krumholz et al. 2013; Li et al. 2013, 2014a). Re¬ 
cent advances in both observational and theoretical studies give 
an opportunity to test the star formation process at small-scales 
(< 1000 AU). 


It has been known for a long time that the presence of mag¬ 
netic fields can drastically change the flow dynamics around low- 
mass stars (e.g., Galli & Shu 1993) and potentially suppresses 
disk formation (e.g., Galli et al. 2006). The latter is due to catas¬ 
trophic magnetic braking where essentially all of the angular 
momentum of the accreting material is removed by twisted field 
lines. Recently, Li et al. (2011) investigate the collapse and disk 
formation from a uniform cloud while Joos et al. (2012) and 
Machida & Matsumoto (2011) performed simulations starting 
with a steep density profile and a Bonnor-Ebert sphere, respec¬ 
tively. They found that rotationally supported disks (RSDs) do 
not form out of uniform and non-uniform cores under strong 
magnetic fields unless the field is misaligned with respect to 
the rotation axis (Hennebelle & Ciardi 2009). Turbulence has 
also been shown to help with disk formation (Santos-Lima et al. 
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2012; Seifried et al. 2012; Myers et al. 2013; Joos et al. 2013; Li 
et al. 2014b). 

In spite of a number of disk formation and evolution simu¬ 
lations, only a few observables have been presented so far. The 
expected observables in the continuum (spectral energy distribu¬ 
tion or SED) from ID and 2D disk formation models have been 
presented in Young & Evans (2005) and Dunham et al. (2010). 
Continuum observables out of 2D hydrodynamics simulations 
with a thin disk approximation have been shown by Dunham 
& Vorobyov (2012) and Vorobyov et al. (2013). However, only 
a handful of synthetic observables from 3D magnetohydrody¬ 
namics (MHD) simulations have been presented in the literature 
(e.g., Commergon et al. 2012a,b). 

Continuum observations probe the dust thermal emission and 
the dust structure around the protostar. However, high spatial 
and spectral resolution molecular line observations are needed to 
probe the kinematical structure as the disk forms. The aim of this 
paper is to present high-spatial (down to 0.1" or 14 AU at a typ¬ 
ical distance of 140 pc) synthetic observations of continuum and 
molecular lines from two of the 3D MHD collapse simulations 
presented in Li et al. (2013). The two simulations differ in the 
initial magnetic field direction with respect to the rotation axis: 
aligned and strongly misaligned where the magnetic field vec¬ 
tor is perpendicular to the rotation axis. The two cases represent 
the two extremes of the field orientation. The synthetic obser¬ 
vations will be compared with those from a 2D semi-analytical 
disk formation model presented in Visser et al. (2009) to inves¬ 
tigate whether the predicted observables differ. The 2D models 
allow us to simplify the different input parameters of the MHD 
simulations into two parameters: sound speed (c^) and rotation 
rate. Rotational transitions of CO are simulated to trace the ob¬ 
servable kinematical signatures. 

Another motivation in simulating CO molecular lines is the 
availability of high-quality spectrally and spatially resolved ob¬ 
servational data toward embedded young stellar objects (YSOs) 
on larger scales (> 1000 AU). Spectrally resolved lines have 
been obtained for low-excitation transitions Ju < 1 = 155 

K) using ground-based facilities (e.g., Jprgensen et al. 2002; van 
Kempen et al. 2009a,b) and higher excited lines up to /u = 16 
(£’u = 660 K) using Herschel-H\F\ (de Graauw et al. 2010) 
in beams of 9-20" (Yildiz et al. 2010, 2013; Kristensen et al. 
2013). Interestingly, San Jose-Garcia et al. (2013) found that 
the C^^O 9-8 lines are broader than the 3-2 lines for low-mass 
YSOs. The observed line widths are larger than that expected 
from the thermal broadening, which indicates a significant con¬ 
tribution of microscopic turbulence or some other forms of mo¬ 
tion, such as rotation and infall. Through the characterization of 
spectrally resolved molecular lines on such physical scales, we 
aim to test the kinematics predicted in various star and disk for¬ 
mation models. 

Another key test of star formation models is to compare the 
predicted mass evolution from the envelope to the star with ob¬ 
servations (e.g., Jprgensen et al. 2009; Li et al. 2014a). Inferring 
these properties toward embedded YSOs is not straightforward 
due to the confusion between disk and envelope. The mass evo¬ 
lution of the disk and envelope can be deduced from millimeter 
surveys combining both aperture synthesis and single dish ob¬ 
servations (Keene & Masson 1990; Terebey et al. 1993; Looney 
et al. 2003; Jprgensen et al. 2009; Enoch et al. 2011). How¬ 
ever, precise determination of stellar masses requires spatially 
and spectrally resolved molecular line observations of the veloc¬ 
ity gradient in the inner regions of embedded YSOs (Sargent & 
Beckwith 1987; Ohashi et al. 1997; Brinch et al. 2007; Lommen 
et al. 2008; Jprgensen et al. 2009; Takakuwa et al. 2012; Yen 


et al. 2013). Here we apply a similar analysis on the synthetic 
continuum and molecular line data as performed on the observa¬ 
tions to test the reliability of the inferred masses. 

This paper is structured as follows. Section 2 describes the 
simulations and the radiative transfer method that are used. The 
synthetic continuum images are presented in Section 3. Section 4 
presents the synthetic CO moment maps and line profiles for the 
different simulations. The results are then discussed in Section 5 
and summarized in Section 6. 

2. Numerical simulations and radiative transfer 

2.1. Magneto-hydrodynamical simulations 

We utilize the 3D MHD simulations of the collapse of a 1 M© 
uniform, spherical envelope as described in Li et al. (2013). The 
envelope initially has a density of po = 4.77 x 10“^^ g cm“^, a 
solid-body rotation of Qq = 10“^^ Hz, and a relatively weak, 
uniform magnetic field of = H (T = 10 where A is the 
dimensionless mass-to-fiux ratio). The details of the simulations 
can be found in Li et al. (2013). 

A snapshot of two simulations dAt = 3.9 x 10^^ s = 1.24 x 
10^ years is used. This corresponds to near the end of the Stage 0 
phase of star formation where almost one half of the initial core 
mass has collapsed onto the star (Robitaille et al. 2006; Dunham 
et al. 2014). The difference between the two simulations is the 
tilt angle between the rotation axis and the direction of initial 
magnetic field vector, Oq. One simulation starts with an initial tilt 
angle Oq = 0° in which a pseudo-disk forms but not an RSD. The 
other simulation starts with an initial tilt angle of Oq = 90° in 
which an RSD forms (see Eig. 1 and Eigure 1 in Li et al. 2013). 

The RSD simulation (left) forms a flattened structure with 
number gas densities > 10^-^ cm“^ in the inner 300 AU ra¬ 
dius. In the region r > 100 AU, the magnitude of the radial and 
azimuthal velocities are within a factor of 2 of each other. The 
radial velocities nearly vanish in the inner 70 AU radius (see 
Eig. 3). The streamlines in the RSD simulation show a coher¬ 
ent flattened rotating component (see Eig. 2). In the case of the 
pseudo-disk simulation, number densities of > 10^-^ cm“^ 
encompass r < 700 AU regions, which is a factor of 2 larger than 
the RSD simulation. An outflow cavity is present in this simu¬ 
lation with an expanding velocity field as shown in Eig. 1 (Left) 
of Li et al. (2013). The cavity is more evacuated compared with 
that in the simulation simulation that forms a rotationally sup¬ 
ported disk although still not a canonical definition of a cavity. 
The magnitude of the radial velocities are much larger than the 
azimuthal velocities in the inner r < 300 AU along most 0 direc¬ 
tions. In this simulation, the streamlines show infalling material 
straight from the large-scale envelope onto the forming star. Us¬ 
ing these two simulations with very different outcomes, we can 
investigate the similarities and differences in both continuum and 
molecular line profiles for pseudo-disk and RSD formation in 
3D. 

2.2. Semi-analytical model 

Eor comparison, synthetic images from 2D semi-analytical ax- 
isymmetric models of collapsing rotating envelope and disk for¬ 
mation as described in Visser et al. (2009) with modifications 
introduced in Visser & Dullemond (2010) and Harsono et al. 
(2013) are also simulated. These models are based on the col¬ 
lapse and disk formation solutions of Terebey et al. (1984), and 
Cassen & Moosman (1981) including a prescription of an out¬ 
flow cavity. The disk evolution follows the tr-disk formalism 
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Table 1. Stellar {M^), envelope (Menv), and disk (M^) masses for the three simulations. is the extent of the rotationally supported disk for each 
simulation. 


Model 

M* 

[Mq] 

-^env 

[Mo] 

Q 

[s-i] 

A 

Md 

[Mol 

[AU] 

3D MHD RSD 

0.38 

0.29 

10“*^ 

10 

0.06 

250-300 

3D MHD No RSD 

0.24 

0.35 

io-‘3 

10 

0.13* 


2D RSD 

0.35 

0.32 

10-'^ 


0.04 

65 


Notes. Mass of pseudo-disk is the sum of the regions with number densities nm > 10^ ^ cm no radius is tabulated for this case. 
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Fig. 1. Density and dust temperature structures in the inner 200 AU radius for all three simulations at t ~ 1.2 x 10^ years. Temperature contours 
at 25 K and 40 K are indicated in the right panels. Top: A vertical slice (R - z slice at 0 = 0 where R denotes the cylindrical radial coordinate) 
of the 3D MHD simulation of RSD formation. Middle: A vertical slice of the 3D MHD simulation of a pseudo-disk (No RSD). Bottom: 2D 
semi-analytical disk formation model. The red lines highlight the region of the stable RSD. 


as described in Shakura & Sunyaev (1973) and Lynden-Bell & 
Pringle (1974). The disk surface is defined by hydrostatic equi¬ 
librium as described in Visser & Dullemond (2010) and is as¬ 
sumed to be in Keplerian rotation. In order to compare with the 
MHD simulations, we consider the collapse of 1 M©, = 0.26 

km s“\ and Qq = 10“^^ Hz core. The sound speed in this case is 
higher than that used in Li et al. (2013), which affects the final 
disk properties at the end of the formation process. The synthetic 
observables are produced at ^ = 3.9 x 10^^ s. The bottom of Fig. 1 


shows the physical structure of the 2D semi-analytical model at 
the time when a ~ 65 AU radius RSD is present. 

A major difference between the 2D semi-analytical axisym- 
metric model and the 3D simulations is the outflow cavity. The 
photon propagation is still treated in 3D. Although outflowing 
gas is present in the pseudo-disk simulation (No RSD), the cav¬ 
ity remains filled with high number density (10^ cm“^) gas while 
lower density (10^“^ cm“^) gas occupies the cavity in the 2D 
model. The outflowing gas generated from angular momentum 
conserving gas in the pseudo-disk model has relatively low ve- 
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Fig. 2. Velocity streamlines for the two MHD simulations in Li et al. (2013). The simulations are rendered with paraview (http://www. 
paraview.org). Left: 3D MHD simulation of a collapsing uniform sphere with the magnetic field vector perpendicular to the rotation axis. 
Right: The simulation with magnetic field vector aligned with the rotation axis. The color of the arrows indicate the azimuthal velocity vector, 
in km s“^ The solid black arrows indicate the general stream lines in the two simulations. The dark shaded colors indicate the density isocontours 
of nH 2 = 10^-^ cm“^ up to 10^^-^ cm“^ (red). 



Fig. 3. Radial {solid) and azimuthal velocities {dashed) at the midplane 
{6 - n 12, (p - Q) for the three different models: 3D RSD MHD simula¬ 
tion {blue), 3D No RSD {red), and 2D model {green). 

locities such that it does not clear out the cavity. As a result, the 
dust temperature along the cavity wall is higher in the 2D model 
due to the direct illumination of the central star. This is readily 
seen in the 40 K contour in Fig. 1 where it is elongated in the 
z direction in the 2D case. However, we show in Section 4 that 
this difference does not affect the results of this paper. 

2.3. Rotation ally supported disk sizes and masses 

The extent of the RSD in the 2D semi-analytical model is defined 
by hydrostatic equilibrium. Using these properties, the RSD is 
defined by a region with densities > 10^-^ cm“^ and azimuthal 
velocities V(p > 1.2 km s“^. The disk evolution in 2D follows 
the alpha-disk formalism with a = 10“^, which in turns define 
the radial velocities as VrlVcj) ~ a{HIR)^ ~ 10“^ where H is the 
disk’s scale height. Similar criteria are used to extract the extent 
of the RSD in the 3D simulation with the additional constraint 
of Note that we apply the criteria from the 2D model to 

define the RSD in the 3D simulation, thus it is is not necessar¬ 


ily in hydrostatic equilibrium. This can be clearly seen in Fig. 3 
for the 3D simulation where Vrive does not satisfy the hydro¬ 
static disk criterion along a few azimuthal angles. Using these 
criteria, an RSD up to 260 AU is found in the mis-aligned simu¬ 
lation. The extent of the 3D RSD is 300 AU if > 1 km s“^ is 
used. With the former criterion, the disk masses contained within 
such a region are 0.06 M© for the 3D RSD and 0.04 M© for the 
2D RSD surrounding 0.38 M© and 0.35 M© stars, respectively 
(see Table 1). As Fig. 3 shows, the radial velocity component 
of the pseudodisk is a significant fraction of the Keplerian ve¬ 
locity. Due to such high radial velocities, the surface density of 
the pseudodisk remains low. However, the total mass of the high 
density regions (^h 2 > 10^’^ cm“^) is a factor of 2 higher than 
the RSD mass. 


2.4. Observables and radiative transfer 

The first step before producing observables is the calculation of 
the dust temperature structure, which is critical for the molecular 
abundances since it controls the freeze-out from the gas onto the 
dust. The dust temperature is computed using the 3D continuum 
radiative transfer code RADMC3D^ with a central temperature 
of 5000 K, which is the typical central temperature in the 2D 
semi-analytical models at around the end of Stage 0 phase. The 
central luminosity is fixed at 3.5 L© for all models. The dust 
opacities used are those corresponding to a mix of silicates and 
graphite grains covered by ice mantles (Crapsi et al. 2008). The 
3D MHD simulation from Li et al. (2013) has an inner radius of 
10^"^ cm (6.7 AU) while the semi-analytical model has an inner 
radius of 0.1 AU. The simulated observables presented here are 
not sensitive to the physical and chemical structures in the inner 
10 AU radius. Thus, these differences do not affect the conclu¬ 
sions of this paper. The same opacities and central temperature 
are adopted for all simulations in order to focus on the general 
features of the observables. The gas temperatures are assumed 
to be equal to the dust temperatures, which is valid for the op- 


^ http://www.ita.uni-heidelberg.de/~dullemond/ 
software/radmc-3d 
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tically thin lines simulated here that trace the bulk mass where 
T^gas ~ T^dust (Doty et al. 2002; Doty et al. 2004). 

CO abundance. In this paper, we concentrate on simulating 
CO molecular lines of / = 2-1, 3-2, 6-5, and 9-8. For sim¬ 
plicity, the abundance is set to a constant value of 10“"^ 
with respect to H 2 except in regions with Tdust < 25 K, where 
it is reduced by a factor of 20 to mimic freeze-out (Jprgensen 
et al. 2005; Yildiz et al. 2013). We adopt constant isotopic ra- 
tios of '2C/'3C= 70, *®0/'^0= 540, and '^0/'’0= 3.6 (Wilson 
& Rood 1994) to compute the abundance structures of the iso- 
topologs. 

Synthetic images. This paper presents synthetic continuum 
maps at 450, 850, 1100, and 1300 yum. The images are rendered 
using RADMC3D with an image size of 8000 AU at scales of 
5 AU pixels. They are placed at a distance of 140 pc. Synthetic 
images at inclinations of 0° (face-on; down the z-axis), 45°, and 
90° are produced. The latter option is included because one of 
the claimed embedded disk sources is close to edge-on (/ ~ 90°, 
L1527 in Tobin et al. 2012). For the synthetic molecular lines, 
the local thermal equilibrium (LTE) population levels are com¬ 
puted using the partition functions adopted from the HITRAN 
database (Rothman et al. 2009). LTE is a good assumption be¬ 
cause the densities in the simulations are greater than the critical 
densities of the simulated transitions. Non-LTE effects may play 
a minor role for synthetic /u ^ 6 lines from the pseudo-disk sim¬ 
ulation due to its lower densities relative to the other two simu¬ 
lations. The line optical depth (r^) is also not expected to play 
a role since the focus of this paper is on the minor isotopologs 
and on the kinematics dominated by the line wings where their 
line optical depths are lower than at the line center. The proper¬ 
ties of the molecules and A^\) are taken from the LAMDA 
database (Schoier et al. 2005). Only thermal broadening is in¬ 
cluded in simulating the molecular lines without any additional 
microturbulence. The image cubes are rendered at a spectral res¬ 
olution of 0.1 km s“^ covering velocities from -7.5 to 7.5 km 
s“^ In order to simulate observations, the synthetic images are 
then convolved with Gaussian beams between 0.1" to 20". The 
convolution is performed in the Fourier space with normalized 
Gaussian images. 

3. Continuum 

3.1. Images and prospects for ALMA 

Continuum images are rendered at four wavelengths and viewed 
at three different inclinations. Figure 4 presents the synthetic 
450, 850, 1100 and 1300 pm continuum images at an inclina¬ 
tion of 45° for the two 3D MHD simulations. The left panels 
present the images from 3D RSD simulation while the right pan¬ 
els show images from the pseudo-disk simulation (labeled as No 
RSD). The features produced during the collapse are clearly visi¬ 
ble in the 450 pm map; most of them are two orders of magnitude 
fainter at 1300 pm. The spiral structure in the 450 pm image is 
due to magnetically channelled, supersonically collapsing ma¬ 
terial on its way to the RSD, rather than a feature of the RSD 
itself. 

One of the aims of this paper is to investigate whether 
these features are observable with current observational lim¬ 
its and what is possible with future Atacama Large Millime¬ 
ter/submillimeter Array (ALMA) data. With the full ALMA, a 
sensitivity of (T=0.5 mJy bm“^ (bm = beam) at ~ 450 pm and 


RSD 


No RSD 



-2 0 2 -2 0 
ARA n 


Fig. 4. Continuum intensity maps of 450, 850, 1100, and 1300 pm for 
the MHD disk formation simulation (left panels) and pseudo-disk (No 
RSD, right panels) at i = 45° at 5 AU pixels. Note the large dynamic 
range needed to see all of the structures. 


0.05 mJy bm“^ at 1100 pm can be achieved at spatial resolutions 
< 0.1" for 30 minutes of integration (1.8 GHz bandwidth). In ad¬ 
dition, a high dynamic range of > 1000 (cr = 0.001 x 5 peak) can 
also be achieved. Figure 5 presents the images convolved with 
a 0.1" beam for the two 3D MHD simulations (right-most pan¬ 
els). The color scale indicates the full range of emission, while 
the red solid lines show the region above 3cr where cr is either 
dynamically limited to 1000 or to the cr values listed above. This 
simply means that the solid red lines indicate the detectable fea¬ 
tures. With a combination of high dynamic range and sensitivity, 
the features of the collapse are observable and distinguishable 
in both 450 pm and 1100 pm continuum maps. At a spatial res¬ 
olution of 0.1", most of the emission at 1100 pm is due to the 
rotationally supported disk with a small contribution from the 
surrounding envelope. 

As Figs. 4 and 5 show, both the strength of the features and 
their extent change with wavelength, as indicated by the red line 
contours. With a resolution of 0.1", the extent of the detectable 
emission decreases from 1.5" at 450 pm to < 1" at 1100 pm for 
the RSD case. At long wavelengths, the dust emission is given 
by ly ~ X (1 “ Since the dust emission is optically 

thin at long wavelengths, the intensity is cx Ldust^'^^dust- The vari¬ 
ables that depend on position are Ldust and Tdust ^ Pdust^v The 
frequency dependence of the opacity is Ky oc for the adopted 
opacity table. The extent of the detectable emission depends on 
these quantities. At one particular position, the ratio of the emis¬ 
sion at the two wavelengths is simply Asoyum/Aiooyum Thus, 

the predicted difference in size at the two wavelengths is due to 
the frequency dependence of the emission. 
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Fig. 5. Convolved continuum maps in the inner 5" of 450 (left) and 1100 (right) pm at i = 45°. The images are convolved with 1", 0.5", and 
0.1" beams as indicated in each panel. For each panel, the top row shows synthetic image of RSD simulation and the pseudo-disk simulation is 
shown in the bottom row. The color scale presents the full range of the emission above 10"^ mJy/bm, not all of which may be detectable. Solid red 
contours are drawn from 3cr up to maximum at 6 logarithmic steps where 1 cr is 0.01 X the maximum (dynamical range of 100) with a minimum 
at 0.5 mJy/bm for the 1" and 0.5" images. For the images convolved with a 0.1" beam, the red contours are drawn with a minimum noise level of 
0.5 mJy/bm for 450 pm and 0.05 mJy/bm at longer wavelengths with a dynamic range of 1000, as appropriate for the full ALMA. 


For ALMA early science observations (cycles 0 and I), the 
capabilities provided a dynamic range only up to 100 at a spatial 
resolution of ~ 0.5". Figure 5 presents synthetic 450 and 1100 
pm images convolved with 1" and 0.5" beams, compared with 
the 0.1" beam images. The color scale again indicates the full 
range of emission, while the red solid lines now show the re¬ 
gion above 3cr where the noise level, cr, is dynamically limited 
to 100 with a minimum of 0.5 mJy/beam at both wavelengths. 
The red lines again present the observable emission. The 450 
pm images convolved with a 1" beam show an elongated flat¬ 
tened structure for both RSD and pseudo-disk simulations and 
are therefore indistinguishable. Most of the emission at 1100 pm 
and longer wavelengths is not detectable at the assumed noise 
level of 0.5 mJy bm“^ A similar result is found after convolu¬ 
tion with a 0.5" beam. Although the synthetic observations from 
the pseudodisk indicate a ‘cometary’ structure, this may be af¬ 
fected by the presence of outflow cavity which is absent in the 
case of RSD. Thus, the full ALMA capabilities are needed to 
distinguish models based on continuum data only. 

3.2. Inclination effects 

The images for a face-on system (i = 0°) are similar to those 
of the moderately inclined system (i = 45°) presented in Fig. 5. 
For low inclinations between 0 to 45°, the continuum images 
only change slightly in terms of absolute flux density and show 
a small elongation due to orientation. 

In contrast, the two simulations rendered at edge-on (i ~ 90°) 
geometry exhibit similar compact components even at the high¬ 
est angular resolution (Fig. 6). They both show an elongated flat¬ 
tened disk-like emission similar to the 2D semi-analytical model. 
This signature suggests an RSD, but it is due to the pseudo-disk 
produced in the 3D MHD simulation whose initial magnetic fleld 
vector is aligned with rotation axis. The peak continuum emis¬ 
sion of the pseudo-disk is a factor of 10 lower than that of the 
other two models while it is similar between the 3D RSD and 
the 2D semi-analytical models (difference of < 10%). The main 


S850^xlxl [Jy/bm] 

-4 -4 -3 -2 



ARA [”] 


Fig. 6. Synthetic 850 pm continuum images convolved with 0.5" (top) 
and 0.1" (bottom) beams for the three simulations viewed at 90° (edge- 
on): RSD formation (left), pseudo-disk (center), and 2D semi-analytical 
model (right). The red solid lines are drawn at the same contours as in 
Fig. 5. 

difference is the extent of the elongated emission where the 2D 
model predicts a very compact (~ 1" radius) structure while the 
pseudo-disk component shows an extended flattened structure 
(< 2"). This illustrates the difflculties in testing disk formation 
models for highly inclined systems based solely on continuum 
data. 

4. Molecular lines 

The continuum emission arises from thermal dust emission and 
does not contain kinematical information. As shown in the pre¬ 
vious section, a compact flattened structure is expected in the 
continuum maps in the inner regions of all three models. Kine- 
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matical information as contained in spectrally resolved molecu¬ 
lar lines is essential to distinguish the models and to derive stel¬ 
lar masses. The rotational lines of CO isotopologs are used to 
investigate this. 

The rotational lines of ^^CO, C^^O, and are simulated. 
CO is chosen since its abundance is less affected by chemical 
evolution during disk formation. We do not investigate ^^CO 
lines since they are dominated by the entrained outflow material 
and are optically thick. The isotopolog lines are more optically 
thin and are expected to probe the higher density region where 
the disk is forming. These predicted spatially resolved molec¬ 
ular line maps can be compared with ALMA data. Moreover, 
high quality spectrally resolved CO isotopolog lines probing the 
larger-scale envelope toward low-mass embedded YSOs have 
been obtained with single-dish telescopes (see §1). The charac¬ 
terization of these C^^O and C^^O line profiles provide a test for 
the kinematical and density structures of the collapsing proto¬ 
stellar envelope on larger scales (e.g., Hogerheijde et al. 1998; 
Jprgensen et al. 2002). 

4.1. Moment maps: RSDs or not? 


u [km s 

- 2.0 - 1.0 0.0 1.0 2.0 



-10 1-10 1-10 1 
ARA [”] 


be noise limited at the line wings since they are weaker than the 
emission near the line center for typical observations with 1-2 
hours integration time. 

Figure 7 presents synthetic first moment (flux-weighted ve¬ 
locity) maps. The moment zero contours at 20% of peak intensity 
indicate the elongation direction. The flattened density structure 
is oriented in the east-west direction (horizontal) for all three 
models as seen in Fig. 1. Coherent velocity gradients from blue- 
to red-shifted velocity are seen in all three models but not neces¬ 
sarily in the same direction. The presence of an embedded RSD 
is revealed in the 3D RSD and 2D model by the coherent blue- to 
red-shifted velocity gradient in the east-west direction similar to 
the flattened disk structure. On the other hand, in the pseudo-disk 
simulation, the velocity gradient is in the north-south direction 
similar to the continuum image as shown in Fig. 5 (right-most 
panels). Such velocity gradients can be mistaken to be along the 
major axis of the disk without higher spatial resolution and sen¬ 
sitivity data. 

A number of effects conspire to generate a velocity gradient 
along the minor axis of the flattened structure in the pseudo-disk 
simulation. First, since the magnetic braking is efficient in this 
case, the dominant motion of the material is in the radial direc¬ 
tion (see Figs. 2 and 3). Second, the flattened structure in this 
particular simulation has lower-density gas than the RSD sim¬ 
ulation (see Section 2.1). The ^h 2 > 10^'^ cm“^ region extends 
up to 700 AU in size and at an angle with respect to the rotation 
axis (see Fig. 1). Thus, the north-south direction in the pseudo¬ 
disk model shows the infalling material along the streamlines 
connecting the large-scale envelope and the central star. 

At moderate inclinations, the pseudo-disk simulation there¬ 
fore shows a coherent velocity gradient in a more-or-less straight 
north-south fine. The velocity gradient changes to an east-west 
direction at high inclinations. At high inclinations, the observer 
has a direct fine of sight on the high density region shown in 
Fig. 1 and therefore the fine emissions pick up the rotational 
motions of the flattened structure similar to the RSD simula¬ 
tions. Furthermore, the skewness that is present in the moment 
one maps of RSD simulations largely disappears at high incli¬ 
nations for the same reason. Thus, it is difficult to separate the 
envelope from the disk for high inclinations (almost edge on) 
from moment one map alone. 


Fig. 7. Moment one maps of 2-1 {top) and 6-5 {bottom) for the 
three simulations viewed i = 45° convolved with a 0.1" beam. The solid 
line shows the 20% intensity contour of the moment zero {f S ydv) peak 
map to indicate the flattened structure. The black dashed lines indicate 
direction at which the PV slices are constructed and the direction of the 
elongation in the moment zero maps. 

Observationally, the kinematical information of the infalling 
envelope and RSD is inferred through moment maps. Elongated 
moment zero (velocity integrated intensity) maps give an indica¬ 
tion of the presence of a flattened structure that is associated with 
a disk. Meanwhile, coherent velocity gradients in the moment 
one (velocity weighted intensity) map may point to a rotating 
component. Analysis of synthetic moment maps of the simula¬ 
tions are presented and compared in this section. We focus on 
presenting the synthetic ‘interferometric’ maps of optically thin 
and fines by convolving the image cubes with a 0.1" 
Gaussian beam. The construction of moment maps only takes 
into account emission > 1% of the peak emission. This translates 
to a noise level of 0.3% of the peak emission. Although ALMA 
can achieve a dynamic range of > 500, it is more likely that the 
molecular fines of minor isotopologs from low-mass YSOs will 


4.2. Velocity profiles 
4.2.1. PV cuts 

Observationally, the presence of embedded Keplerian disks is of¬ 
ten established by constructing position velocity (PV) diagrams 
along the major axis of the system as seen in moment zero maps. 
In theory, the PV analysis is straight-forward. It is symmetric in 
both position and velocity space (4 quadrants are occupied) if 
the system is infall dominated (e.g., Ohashi et al. 1997; Brinch 
et al. 2008). The symmetry is broken if rotation is present and the 
emission peaks are shifted to larger offsets corresponding to the 
strength of the rotational velocities (2 quadrants are occupied). 

Figure 8 presents synthetic PV diagrams along the major axis 
of the disk where it corresponds to the direction of the blue- to 
red-shifted velocity gradient. For the images in Fig. 7, an east- 
west slice (horizontal) is taken for the RSD simulations, while 
a north-south (vertical) slice is adopted for the case of the sim¬ 
ulation without an RSD (No RSD). These slices are not exactly 
the major axis of the moment zero map, however these direc¬ 
tions pick up most of the velocity gradient present in the inner 
1". Both and lines are simulated. Interestingly, the 
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RSD No RSD 2D Model 



RSD No RSD 2D Model 



Fig. 8. PV maps of (leff) and {right) along the velocity gradients as seen in Fig. 7 for an inclination of 45° and 0.1" beam. The top panels 
show the J - 2-1 line and the bottom panels show the 6-5 transition. The blue rectangle indicates the radii of RSDs. The dashed lines indicate the 
inclination corrected Keplerian curves associated to the stellar mass. The red color scales show emission from 10% to the peak intensity. 


PV slices suggest that rotational motions are present regardless 
of whether an RSD is present or not. This is most readily seen in 
the PV maps (right of Fig. 8) in which only 2 of the 4 quad¬ 
rants are occupied by molecular emissions for all three models 
at i = 45°. This shows that emission readily picks up the 
rotational motion at small-scales but also that infalling motion 
can be confused with rotation if the wrong direction for the PV 
cut is chosen (No RSD model). 

The PV maps indicate contributions from the infalling 
envelope since the maps are more symmetric than those in C^^O. 
The PV slices of the pseudo-disk simulations indicate an 
infall-dominated structure in which the 4 quadrants are occu¬ 
pied. On the other hand, there is a clear indication of a rotating 
component for the two RSD simulations in which only 2 of the 4 
quadrants are filled at small radii. This suggests that spatially and 
spectrally resolved lines can distinguish between a pseudo¬ 
disk and an RSD. 

Figures 8 compares the PV maps of the 2-1 and 6-5 transi¬ 
tions. Most of the emission in the 6-5 transition occupies only 
2 of the 4 quadrants indicating signatures of rotational motions, 
whereas the 2-1 lines also show some emission in the other 2 
quadrants from larger scales. This is a clear indication that the 
6-5 line is a better probe of the rotational motions in the in¬ 
ner 100 AU than the 2-1 line. Yet, a combination of spatially 
and spectrally resolved molecular lines are needed to confirm 
the presence of embedded rotationally supported disk. 

4.2.2. Peak position-velocity diagrams 

While it is clear that there is indeed a rotating component for 
some models, the question is whether the extent of the Keplerian 
structure can be extracted from such an analysis. A Keplerian 
rotating flattened structure exhibits a velocity profile v oc 
where r is the distance from the central source. These positions 
are either determined from fitting interferometric visibilities of 
each velocity channel (Lommen et al. 2008; Jprgensen et al. 
2009) or determination of the peak positions in the image space 
(Tobin et al. 2012; Yen et al. 2013). We here determine the peak 
positions directly in the image space to assess whether a velocity 
profile is visible in the synthetic molecular lines. 


Peak positions are determined for each of the velocity chan¬ 
nel maps for each molecular line for the red and blue-shifted 
components separately taking into account channels whose peak 
fiux density {Sy in Jy bm“^) are > 1% of 5max- They are sub¬ 
sequently rotated according to the direction of the velocity gra¬ 
dient. If an RSD is present, the peak positions of both red- and 
blue-shifted velocities are expected to follow the Keplerian ve¬ 
locity profile {v oc The combination of infalling rotating 

envelope and RSD, which exhibits a skewness in the moment 
one map, is expected to show a steeper velocity profile {v oc r~^) 
(Yen et al. 2013). However, at high velocities, the peak positions 
of the red- and blue-shifted velocities can be misaligned at scales 
of 5 AU due to the limited spatial resolution. The disk radius is 
determined by minimizing the difference (~ 10%) between the 
best-fit stellar mass inside and at the disk radius. 

The peak position-velocity diagrams (PPVs) for 6-5 
and 2-1 are shown in Fig. 9. These lines are chosen be¬ 
cause they represent two observational extremes in terms of ex¬ 
citation and optical depth. There is a clear distinction between 
the RSD simulations and a pseudo-disk. The velocity profile of 
the pseudo-disk is much flatter than that expected from an RSD. 
Thus, spatially and spectrally resolved molecular lines observa¬ 
tions can clearly differentiate between an RSD and a pseudo¬ 
disk. 

All three PPVs indicate a velocity profile close to i; cx 
(for the pseudo-disk see 2-1 PPV in the inner 20 AU). This 
reflects the fact that the inner 300 AU of the models is domi¬ 
nated by velocity structure that is proportional to which is 

both the radial velocity ^r'mfaii ^ ^2^^* j and angular velocity 

^^rot ^ (Brinch et al. 2008). From such characterization, 

the stellar masses can be calculated and indicated in the top right 
corner of Fig. 9. In general, the best-fit stellar masses in the case 
of RSD are within 30% of the true stellar masses tabulated in 
Table 1. This is not so for the 2D model in the 6-5 and 
also 6-5 (not shown) due to the fact that the inner flattened 
envelope is warm (> 40 K) and dense. 

Another parameter that one would like to extract is the disk 
radius, as indicated by the vertical solid line in Fig. 9. The 
break at is readily seen in the 2D model at ~ 40 AU for both 
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Fig. 9. Peak position-velocity diagrams of 
6-5 {top) and 2-1 {bottom) at 
i = 45° after convolution with a 0.1" beam. 
The red and blue symbols correspond to the 
red- and blue-shifted velocity components, 
respectively. The different panels present the 
PPV for the different simulations. Vertical 
solid lines show the disk radii extracted for 
the three models and dashed lines show 
the steep velocity profile {v oc r~^) while 
the solid cyan lines indicate the Keplerian 
curves. The stellar masses are indicated in 
the top right of each panel. The offset be¬ 
tween blue- and red-shifted points are due 
to the limited spatial and spectral resolution 
at the high velocities. 


6-5 and 2-1. For the case of the 3D MHD simulation 
(RSD), the best-fit radius varies between 100 and 300 AU. It 
also exhibits a steep velocity profile {v oc r at radii >100 
AU. The large range of disk radii is due to the envelope emission 
overwhelms the molecular emission because CO is frozen out in 
the cold part of the disk at those large radii. The issue of disk 
versus envelope emission becomes apparent in the case of large 
embedded disk {R^ >100 AU). 

The comparison shows that there is a clear distinct PPV pro¬ 
file associated with RSD formation from to r~^. A steep 
velocity profile {v oc r“^) is absent in the pseudo-disk simula¬ 
tion. This seems to indicate that such a steep velocity profile 
describes on-going RSD formation based on the given simula¬ 
tions. A pseudo-disk is characterized by a flat velocity profile 
in the inner regions. Furthermore, the PPV method can simulta¬ 
neously derive the stellar mass and the extent of the RSD while 
separating the infalling rotating envelope from it. With respect to 
differentiating between RSD and non-RSD, PPV is a better tool 
than PV-diagrams. 

4.3. Single-dish line profiles 

The previous sections focus on features at small-scales as ex¬ 
pected from interferometric observations. The next assessment 
is to compare the synthetic molecular lines with single-dish ob¬ 
servations which probe the physical structure of the large-scale 
envelope on scales up to a few thousand AU. The image cubes 
are convolved with 3 different beams: 9", 15", and 20". These 
different beams are typical for single-dish CO observations using 
the JCMT (15"), Atacama Pathfinder Experiment (APEX, 9"), 
and Herschel (20"). Figure 10 presents the synthetic CO lines 
(/u = 3, 6, and 9) for the two MHD simulations viewed face- 
on {i ~ 0°) convolved with a 9" beam. The face-on orientation 
is considered first to compare with the line profiles in Harsono 
et al. (2013) for the 2D simulation. The low-lying transitions {J^ 
= 3) probe the kinematics in the large-scale envelope. 

Double peaked line profiles are present in the ^^CO and C^^O 
3-2 regardless whether an RSD is present or not. For the ^^CO 
line, an inverse P-Cygni line profile is seen due to the coherent 


infalling material onto the disk while a P-Cygni profile is associ¬ 
ated to the pseudo-disk, which is tracing the expanding material 
due to outflowing material present in the pseudo-disk simulation. 
Self-absorption causes the double peak in the ^^CO 3-2 line due 
to optical depth (typically, tl > 5) at line center whereas it is 
weakly affecting the C^^O 3-2 line. 

It is interesting to note that there is no significant difference 
in the 6-5 lines between a simulation that forms an RSD versus a 
pseudo-disk. This transition {E^ -110 K) traces the dense warm 
gas where a large fraction of the emission comes from > 40 K 
gas (Yildiz et al. 2010). The P-Cygni line profile is still visible 
in the ^^CO line, however it is not significant in the C^^O line. 
The lines are also not Gaussian with significant wing emission 
extending up to +2 km s~^. 

The line profiles for the semi-analytical models within a 9" 
beam are shown in the bottom of Fig. 10. They are signifi¬ 
cantly different from the 3D MHD models, which arises from 
the prescribed velocity structure. In Harsono et al. (2013), an 
additional microturbulent broadening of 0.8 km s“^ was added, 
which results in Gaussian line profiles consistent with the ob¬ 
served single-dish CO line profiles. However, in this paper, we 
have not included the additional broadening term in order to in¬ 
vestigate the emission arising from the true kinematical informa¬ 
tion. The peaks are more prominent than those in the 3D MHD 
simulations due to a jump between the velocity structures of the 
RSD component and the infalling envelope. In general, the 2D 
semi-analytical models produce significantly broader lines and 
significant variations between the CO isotopologs and transitions 
compared with 3D simulations because of a warmer disk and 
outflow cavity wall (see Section 2.2) which allow for stronger 
wing emissions. 

4.3.1. Inclination effects 

The simulated lines viewed face-on may not pick up all of the 
dynamics of the system. Figure 11 shows how the line profiles 
change with inclination for the three different simulations within 
a 15" beam. The lines become broader with increasing inclina¬ 
tion as they readily pick up the different velocity components. 
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Fig. 11. line profiles for the different simulations viewed at i - 45° 
and i = 90° in 15" beam. 

Table 2. FWHM (FWZI as defined at 10% of the peak emission) in 
km s~^ of the and lines within a 20" beam for the MHD 
simulations viewed at i - 45°. 


/ 

D 

RSD 

CO 

No RSD 

C 

RSD 

No RSD 

3-2 

0.8 (1.8) 

1.1 (2.2) 

0.9 (1.8) 

0.9 (1.9) 

6-5 

1.1 (2.3) 

1.0 (2.2) 

1.3 (2.3) 

0.9 (2.2) 

9-8 

2.1 (3.7) 

1.2 (3.9) 

2.3 (4.0) 

1.0 (3.9) 


Fig. 10. Top: i^CO, C^^O, and 3-2 {top\ 6-5 {middle), and 9-8 
{bottom) spectra at i ~ 0° within a 9" beam. The blue line shows the 
synthetic line from simulation with an RSD while the red line is the 
simulation without an RSD. Bottom: Spectral lines convolved with a 
9" beam simulated from 2D semi-analytical model viewed at face-on 
orientation. 


the other hand, in the case of pseudo-disk formation, both radial 
and azimuthal velocities are of equal importance. The origin of 
the line broadening therefore depends on whether or not an RSD 
is forming. If an RSD is indeed forming, the 9-8 is broad¬ 
ened by rotational motions at moderate and high inclinations; at 
low inclinations, infall dominates the broadening. 


The /u = 3 lines exhibit inverse P-Cygni line profiles that are 
associated with infalling gas. Meanwhile, the higher J transi¬ 
tions show more structured line profiles compared with the sys¬ 
tems viewed face-on. The 6-5 lines are double-peaked in 
both cases of RSD formation while it is single-peaked for the 
pseudo-disk model. On the other hand, the 9-8 line is signifi¬ 
cantly broader than the low-/ lines reflecting the complexity of 
the dynamics of the warm dense gas. 

4.3.2. Origin of the line broadening 

To investigate the source of the line broadening, a set of molec¬ 
ular lines are simulated with zero radial velocity {v^ = 0 km s“^) 
and another with zero azimuthal velocity = 0 km s“^). For 
the MHD simulation with RSD, the FWHM value of the 
9-8 line decreases to < 0.5 km s"^ at all inclinations without any 
azimuthal velocity component. Such a decrease is not dramatic 
for face-on orientation, however it is more than a factor of 3 for 
intermediate {i ~ 45°) and high inclination {i > 75°) cases. On 


4.3.3. Line widths and comparison with observation 

Molecular line observations are typically characterized by their 
peak flux densities (or intensities), FWHM, and integrated line 
flux densities. While the peak flux densities and integrated line 
fluxes depend on the adopted physical and chemical structure, 
their FWHM should reflect the general kinematics that are 
present in the system. In this paper, we focus on the comparison 
of FWHM and full-width at zero intensity {FWZI) determined 
at 10% of the peak with observations. A 10% cut-off is chosen 
since most single-dish observations do not reach higher signal- 
to-noise, especially for the higher-/ lines. 

These values are calculated for and lines of the 
different models from the convolved image cubes. The FWHM 
and FWZI within a 20" beam are listed in Table 2 at moderate 
inclination {i = 45°) comparing the two 3D MHD simulations. 
Within such a large beam, their values for and are 
similar. The FWHM values do not necessarily increase between 
/u =3 and 6, in contrast with the FWZI values (see Fig. 12). 
This is expected since the wing emissions are much lower than 
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Fig. 12. FWHM values of the observed (gray shaded region) and simu¬ 
lated (symbols) rotational lines of and as indicated in a 20" 
beam. The observed values are taken from J0rgensen et al. (2002), van 
Kempen et al. (2009a), and San Jose-Garcia et al. (2013). The FWHM 
of synthetic lines are indicated at specific inclinations as indicated by 
the symbols: circles for 0°, triangles for 45°, and squares for 90°. The 
comparison of the FWZI values are specifically for the 3-2 line. 
The different colors represent the different type of simulations: 3D RSD 
(RSD), 3D pseudo-disk (No RSD), and 2D semi-analytical model (2D 
Model). 


the peak because the emitting region is much smaller than the 
beam (beam dilution). 

Considering all inclinations (0°, 45°, and 90°), the FWHM 
values of the low-/ CO lines are similar between the 3D simu¬ 
lations and 2D semi-analytical models. On the other hand, the 
line widths of the high-/ lines differ significantly. This is ex¬ 
pected since most of the 2-1 and 3-2 emission originates from 
the large-scale envelope, which is similar in terms of kinemat¬ 
ics in the three simulations. However, the high-/ lines originate 
from the warm inner regions in which the three simulations show 
different velocity structures. 

Figure 12 presents the comparison between the simulated 
and observed lines. The observed line widths toward a sample of 
Class 0 low-mass YSOs are taken from Jprgensen et al. (2002), 
van Kempen et al. (2009a), and San Jose-Garcia et al. (2013, 
FWHM^ for the narrow component). Sources with known con¬ 
fusion in their line profiles from other nearby sources have been 
excluded. It is clear that the observed line widths of the low- 
/ lines (/u = 2 and 3) are significantly greater than the model 
simulations, by a factor of 2. Since these lines probe the large- 
scale quiescent envelope, the discrepancy between the predicted 
and observed line widths suggests that the large-scale envelope 
is turbulent with FWHM ^ 1 km s“^ or a Doppler of ~ 0.4 km 
s“^ i.e., more than what is included in the current simulations 
that could be due to an interaction with fast outflow, at least in 
part. The comparison of the FWZI values also supports this con¬ 
clusion. 

For the dense and warm gas probed by the higher /u > 6 
lines, the predicted line widths are within the observed range for 
the models that do form an RSD. Moreover, the pseudo-disk (No 
RSD) simulations predict much smaller line widths than those 
observed. This line originates from the inner warm parts that are 
rotationally supported. Our treatment of the inner radius does 
not affect this since an evacuated cavity is absent in the case of 
RSD formation. In fact, if an outflow cavity is allowed to form, 
the 3D RSD predicted line widths would shift upward. Thus, the 


line width of this particular line reflects the true kinematics in 
the inner regions. This comparison may suggest that the large 
observed line widths of 9-8 line indicate that the kinemat¬ 
ics of the warm inner envelope are similar to that of a model that 
forms an RSD. Alternatively, a turbulence with FWHM = 2 km 
s“^ would be required in the inner parts if RSDs are absent. 


5. Discussion 

5.1. Variations with viewing angles 

This paper presents synthetic continuum maps and CO iso- 
topolog lines from 3D MHD simulations and 2D semi-analytical 
disk formation models out of collapsing rotating envelopes. The 
aim is to present signatures that can differentiate between em¬ 
bedded rotationally supported disks (RSDs) and a pseudo-disk. 
Thus far, we have analyzed the continuum and synthetic molec¬ 
ular lines with respect to one viewing angle in the azimuthal di¬ 
rection with 0 = 0°. As shown previously by Smith et al. (2012), 
the line profiles may change with different viewing angles since 
the collapse process is not spherically symmetric. 

In order to look at the general trend with viewing angles, syn¬ 
thetic images at 4 different inclinations (/) from 15° to 150° and 
8 azimuthal angles from 8° to 315° are generated. We concen¬ 
trate the analysis on 850 yum, 1100 yum, C^^O, and images. 
Figure 13 presents the continuum and molecular lines at a few 
viewing angles for the two 3D MHD simulations. The 3D RSD 
formation predicts an observable spiral feature at near face-on 
(i ^ 0°) with high dynamic range (1000). As mentioned earlier, 
this is due to the infalling material on its way to the disk. 

In terms of kinematical signatures, the moment one maps of 
the 2-1 line are compared for the two cases. Both the RSD 
and pseudo-disk simulations predict a coherent velocity gradient 
in the inner 300 AU. However, the velocity gradient is more ro¬ 
bust if an RSD is forming. The pseudo-disk simulation shows a 
velocity gradient in the moment one map that is not necessarily 
corresponding to rotation (see Fig. 1 and Section 4.1). Such a 
direction corresponds to the streamlines of the infalling material 
from the large-scale envelope to the central star. At high incli¬ 
nations (i > 75°), this direction of the velocity gradient shifts to 
east-west direction similar to the RSD simulation, because the 
moment one map is dominated by the rotational motions which 
are in the east-west direction. 

To assess the general predictions for the large-scale enve¬ 
lope, 9-8 spectra within a 15" beam are compared. The 
low-/ lines (/u = 2 and 3) exhibit inverse P-Cygni profiles in¬ 
dicating infalling material in both RSD and pseudo-disk simu¬ 
lations. For the high-/ lines, the 3D RSD simulation predicts a 
broader line than the pseudo-disk simulation in most orientations 
consistent with Fig. 12. 

In summary, the results that are presented in previous sec¬ 
tions are robust and do not depend on the viewing angles. A 
pseudo-disk shows more distinct features in the kinematics in 
the moment one maps that are different from an RSD. A coher¬ 
ent blue- to red-shifted velocity gradient in the inner 1 arcsec- 
onds aligned with the major axis is most likely a signature of 
an RSD or on-going RSD formation. Although these lines are 
simulated assuming LTE conditions, non-LTE effects generally 
decrease the strength of the emission but do not alter the results 
of the kinematics derived from the molecular limes. 
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850 nm 




Fig. 13. Comparison of synthetic 850 /jm, 1100 ^um, C'^O 2-1 moment one maps and normalized C'*0 9-8 spectra for the 3D RSD (left) and 
pseudo-disk {right) MHD simulations viewed at a few orientations. The images and image cubes are convolved with a 0.1" beam. The maps sizes 
are 5" across similar to Fig. 5. The color scales and line contours are defined in Figs. 5 and 7. The normalized 9-8 spectra refer to a 15" 
beam on the same scale as Fig. 11. The inclination (/) and azimuthal angle (0, rotation angle of the observer around the z axis) are representative 
only and not exact. 




Fig. 14. Left: Mean deviation of disk masses (gas -i- dust) derived from 1100 pm continuum and molecular lines fluxes from the true disk mass 
extracted from the simulations. The gray bars indicate the full range of values for the different orientations. Right: Mean deviation of stellar 
masses measured from molecular lines from the true stellar mass. The colors indicate the different simulations: 3D MFID simulations in blue, 3D 
MHD pseudo-disk simulation in red, and 2D semi-analytical model in green. The error bars indicate the standard error of the mean. The shaded 
region shows the area within a factor of 2. 


5.2. Masses and sizes 

5.2.1. Analysis with continuum 

To test models of star and disk formation (e.g., Hueso & Guillot 
2005), the flow of mass from envelope to disk to star with time 
needs to be known. In practice, this means that the masses of 
the central (proto-)star and of its envelope-disk system need to 


be inferred from observations and models at different evolution¬ 
ary stages. A pragmatic way to ‘extract’ the envelope and disk 
masses is to compare the continuum fluxes within a large beam 
with those found at ‘small’ scales. Here, 5" area is taken as the 
‘small’ scale corresponding to 50k/I at a distance of 140 pc as 
probed previously with the SMA. Specifically, Jprgensen et al. 
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(2009)^ gives the following formulae: 


S 5" — 5* disk + ^nv ^ S Qny 

kS 15" = 5 * disk + k^env? 


( 1 ) 

( 2 ) 


where /env is the envelope contribution at ‘small’ scales, Ss" is 
the flux within a 5" area, 5 15 " is the flux within a 15" beam, 
5 env is the total envelope flux, and 5 disk is the total disk flux. 
The set of equations above roughly separates the envelope from 
the disk in terms of continuum flux contribution. There are two 
issues that need to be addressed in utilizing the equations above: 
the envelope fraction (/env) and the ‘small’ scale flux. 

The first issue is presented in §4.1 of Jprgensen et al. (2009). 
They calculated the envelope contribution as determined from 
Pdust ^ spherical envelopes. This yields a maximum con¬ 
tribution of 8% at 850 jum within a scale of 5" compared with 
that at 15". The contribution increases with increasing power- 
law exponent (16% for p oc Note that Jprgensen et al. 

(2009) computed the relative contribution between 5 50 kd at 1 . 1 
mm and 15" in 850 pm. In order to compute the contributions 
within the same wavelengths, these values are scaled assuming 
Fy oc In addition, for a given 2D embedded disk model, 
the projected spherical envelope model is one with an increasing 
exponent with inclination (i.e., p oc r~^-^ as i 90°). In other 
words, for an embedded YSO viewed at face-on, the best repre¬ 
sentative spherical model is one with a rather flat density profile. 
As a result, the envelope contribution from an embedded YSO 
on small scales changes with inclination. 

Another issue is the ‘small’ scale flux at which to determine 
the envelope contribution. If the size is too small (i.e., 1"), the 
envelope contribution is naturally negligible and it is easy to ex¬ 
tract the disk flux. However, the disk flux can be underestimated 
within such a small scale. Furthermore, the flux that is measured 
at each pixel is always a combination of the envelope and the 
disk. Thus, it is more intuitive to calculate the envelope’s contri¬ 
bution at scales larger than the disk (~ 5", 700 AU at 140 pc) in 
order to obtain all of the disk’s flux. 

For this purpose, following the method in the literature, the 
envelope and disk masses are determined from comparing the 
fluxes within a 15" beam and a 5" area. This procedure is per¬ 
formed on the images rendered at the 32 different viewing angles 
as described in the previous paragraphs. A typical constant en¬ 
velope fraction of 8 % at 850 pm (2% at 1100 pm) is adopted 
to see how well such a simple estimation based on spherically 
symmetric envelope model can extract the disk properties of the 
2D and 3D simulations. 

Once the disk and envelope fluxes at 850 pm have been 
obtained, the conversion to dust mass follows Jprgensen et al. 
(2009), Eq. (1), which takes into account the fact that there is 
a distribution of temperatures in the envelope set by the lumi¬ 
nosity of the protostar. The envelope mass therefore scales both 
with distance and bolometric luminosity (Lboi). The bolometric 
luminosities are calculated by constructing the spectral energy 
distribution at the same viewing angles as the synthetic images. 
For the disk, a dust temperature of 30 K is adopted to calculate its 
mass as it is done in observations. These values are compared to 
the masses of different components in the simulations (Table 1). 
The disk mass in the simulation is defined by the velocity struc¬ 
ture as described in Section 2.3. The left-over material within the 
computational box is the envelope mass. 

Figure 14 presents the ratio of disk masses derived from syn¬ 
thetic continuum observations to the true masses tabulated in Ta- 

^ Note that the factor (^) is omitted since the analysis is performed 
on the same wavelength. 


ble 1. In the case of the pseudo-disk, the disk mass is taken to be 
the mass of the region with number densities > 10 ^-^ cm“^. 
The masses inferred from the continuum fluxes are within a fac¬ 
tor of 2 of the true disk mass. There is a little difference in the 
disk mass between 1100 pm and 850 pm for both 3D RSD and 
2D semi-analytical model. Envelope masses agree if scaled to 
the same 15" beam. 

The spread in the inferred disk masses is greater in the 2D 
model than in the 3D simulation. The lower end of the spread is 
occupied by simulations at high inclinations (/ > 75°). In such an 
orientation, the continuum optical depth even at 1100 pm is high 
since the disk is viewed edge-on and, consequently, a large part 
of the disk does not contribute to the observable emission. The 
high end of the spread is when the system is viewed almost face- 
on where the inner flattened inner envelope deviates from the 
spherically symmetric assumption. Finally, the equations above 
underestimate the disk mass from the 3D simulations because a 
Tdust = 30 K was assumed to obtain the disk mass whereas the 
temperature within the disk is < 30 K at large radii. A lower 
temperature ~ 10 K would be best for the 3D RSD simulation 
while it is inclination dependent for the 2D model case. 

The results suggest that disk masses can be well estimated 
from the continuum flux even for the pseudodisk. However, there 
is a clear difference in how the disk is defined in all three differ¬ 
ent cases. In the 3D and 2D RSD, the disk mass is defined by 
its velocity structure such that it is indeed rotationally supported 
while the pseudodisk is defined by density. Hence, in the latter 
case, the difference between the disk and the infalling envelope 
is not well defined. 

In summary, we And that the disk masses as inferred 
from continuum emission in the 3D simulations and 2D semi- 
analytical model are within a factor of 2 of the true value. This 
factor-of -2 is due to the viewing angle and the dust temperature 
within the disk is lower than the assumed 30 K. However, we 
And that the pseudo-disk also indicate similar emission at such 
scales (see also, Chiang et al. 2008). In this case, the disk does 
not corresponds to a Keplerian disk. 


5.2.2. Disk radii and stellar masses 

Analysis of the continuum toward a large sample of sources 
yields disk and envelope masses. In order to fully test the evo¬ 
lutionary models, the stellar masses need to be derived from the 
kinematics of molecular lines. In addition, the disk radii can be 
determined from the velocity profiles (see Section 4.2) and can 
also serve as tests for the evolutionary models since disk radius 
is expected to increase with evolutionary state but to also depend 
on the initial angular momentum of the core (see Fig. 13 in Har¬ 
sono et al. 2014). The extent of the RSD is determined from the 
break in the peak-position velocity diagrams (see Fig. 9) for the 
simulated 2-1, 6-5, 2-1 and 6-5 lines. 

For the 2D semi-analytical models, the average radius from the 
4 lines is 64 + 8 AU as derived from the break. For the 3D MHD 
RSD simulation, the average radius is 230 + 100 AU. Thus, the 
inferred drop-off gives a good estimate of the extent of small 
disks, while a larger disk has a larger error bar associated to it. 

Disk masses can also be calculated from the integrated flux 
densities (f Sydv in Jy km s“^) of molecular lines. In this 
method, the integrated flux densities are extracted within the re¬ 
gion defined by the radius in the previous paragraph. The mass 
calculation is given by (Scoville et al. 1986; Momose et al. 1998; 
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Hogerheijde et al. 1998): 


M, 


gas 


5 9 ^ jq 6_ QjTex) _ 

g-uAui exp (-EalkBTsJ [X/H 2 ] 


Tl 

1 - eXp“^L 



Sydv Mq, 


(3) 


where A^\ is the Einstein A coefficient of the transition, is 
the degeneracy of the upper level, is the upper level en- 
ergy, QiT^^) is the partition function at an excitation tempera¬ 
ture Tex = 40 K, yUH 2 is the mean weight of the gas, [X/H 2 ] is the 
abundance of molecule X with respect to H 2 , tl is the fine opti¬ 
cal depth, and is the Boltzmann constant. A higher excitation 
temperature for the gas is adopted with respect to the dust tem¬ 
perature following the rotational temperature derived from the 
observed lines (Yildiz et al. 2013). The derived masses do 
not strongly dependent on the adopted excitation temperature. 
A constant line optical depth of 0.5 for the and 0.3 for 
the lines is used, which characterizes the average optical 
depth over all velocities. The crucial assumption here is the CO 
abundance, [Y/H 2 ], which is taken to be 10“"^ similar to that in 
Section 2.4. Since the CO abundance will be affected by freeze- 
out in the cooler parts of the disk, this assumption provides a 
lower limit to the disk mass. 

Figure 14 presents the masses obtained from the integrated 
fine flux densities. They show a much larger scatter than the 
masses obtained from the dust continuum flux. The symbols in¬ 
dicate the values appropriate for masses derived at moderate in¬ 
clinations. For the case of RSD formation, the low-end of the 
spread is due to near face-on orientation (30° > i > 150°) in 
which the obtained disk sizes are half of the true RSD and, con¬ 
sequently, the disk mass is lower than average. In the 2D semi- 
analytical model case, the physical structure of the inner enve¬ 
lope plays a role. The conversion between integrated flux to mass 
also assumes that the emission at all velocities is dominated by 
the disk, which is only true for emission at the wings (high ve¬ 
locities > 1 km s“^). However, in the 2D semi-analytical model, 
the inner envelope is highly flattened compared to that of the 3D 
MHD simulation such that high density gas (^h 2 >10^ cm“^) 
surrounds the RSD (see Fig. 1). So, the disk does not dominate 
the integrated flux density at this particular time (see Fig. 8 of 
Harsono et al. 2013) and, consequently, the integrated flux den¬ 
sity is not correlated with its mass. Therefore, a good estimate 
of the disk’s molecular mass can only be obtained if the system 
is oriented at moderate inclinations and the RSD is assumed to 
dominate the molecular emission at all velocities. 

The disk masses derived from the 2-1 line are generally 
higher than that from the 6-5 line. Since the RSD size obtained 
from the 6-5 emission is smaller than that from the 2-1 line (e.g., 
200 vs 100 AU), the integrated fine flux density is extracted from 
a smaller region, which results in a much lower mass. The size 
of the RSD is smaller in the 6-5 emission because the emission 
arises from the dense warm regions in the vicinity of the proto¬ 
star and, thus, its mass. The masses obtained from the integrated 
2-1 line are better estimates for the true disk mass. 

The 3D MHD and 2D semi-analytical RSD formation pre¬ 
dict the same behaviour in terms of disk masses in both 
and lines. The difference of the masses obtained in the 2-1 
and the 6-5 lines are also similar for the two isotopologs. The 
2D semi-analytical model predicts a smaller spread than the 3D 
RSD model in the disk masses obtained from the 2-1 line 
since the high density region is more compact in the 2D case 
and, therefore, the 2-1 integrated flux density obtained in 
the 2D model does not vary as much as in the 3D simulation. 


The comparison between the 2D semi-analytical model and 3D 
MHD shows that the reliability of the observables to trace the 
true masses depends on the physical structure of the inner enve¬ 
lope. 

The stellar masses obtained from the PPV method are shown 
in Fig. 14. They are typically within a factor of 2 of the true 
stellar masses. The best-fit values for synthetic molecular lines 
viewed at / < 15° tend to be more than twice the true stellar mass. 
This is due to the difficulties in obtaining the peak positions at 
high velocities since they are most likely to be below the noise 
level. The 2D semi-analytical model indicates better agreement 
with true stellar masses than those of the 3D model because their 
velocity structures are different. The kinematics of the envelope 
is the same at all viewing angles in 2D. However, there are lines 
of sights that pick up significant sub-Keplerian gas infalling onto 
the disk in 3D as shown in Fig. 3. Therefore, stellar masses as 
obtained from spatially and spectrally resolved molecular obser¬ 
vations are good estimates of the true stellar mass within a factor 
of 2, which is smaller than the uncertainty in the inclination. 

From the analysis of the synthetic observations, spatially 
(< 0.1" at 140 pc) and spectrally (< 0.1 km s“^) resolved op¬ 
tically thin molecular fines at two energy levels (e.g., 2-1 and 
6-5) are required to understand the physical structure of the in¬ 
ner envelope. The disk masses obtained from the two lines are 
similar in the pseudodisk case while they can differ by an or¬ 
der of magnitude in RSD case. The simple PPV analysis directly 
from the data can already differentiate between the RSD and a 
pseudodisk. However, sophisticated modelling tools are needed 
to infer the physical structure of the disk. 

6. Summary and conclusions 

We have presented the observables in continuum and molecular 
fines for two 3D MHD simulations of Fi et al. (2013) and 2D 
semi-analytical models of collapse and disk formation. Snap¬ 
shots of two different MHD simulations of a relatively weakly 
magnetized core {Bq = 11 pG) at the same time after the on¬ 
set of collapse are used. One simulation has an initial mag¬ 
netic field axis aligned with the rotation axis in which a rota¬ 
tional supported disk (RSD) does not form, however a pseudo¬ 
disk and outflowing gas are present. The other MHD simula¬ 
tion starts with the magnetic field axis oriented perpendicular 
to the rotation axis, which results in the formation of an RSD 
(see Figs. 1 and 3). These simulations explore the two extremes 
of the magnetic field orientation. The synthetic observables are 
then compared to a 2D semi-analytical model without magnetic 
field (Visser et al. 2009) with similar initial conditions (1 M© and 
Qo = 10“^^ Hz). Accurate dust temperatures are calculated using 
the 3D continuum radiative transfer tool RADMC3D with the 
same dust opacities and central temperature for all three mod¬ 
els. Continuum images and thermalized CO molecular lines are 
produced using the same radiative transfer code and method. 
Freeze-out of CO onto dust grains is included. This paper fo¬ 
cuses on presenting similarities and differences in the predicted 
observables. The main results and conclusions are as follows. 

- Synthetic continuum images of the two MHD simulations 
and 2D semi-analytical model indicate that a spatial resolu¬ 
tion of 14 AU and high dynamic range (1000) are required to 
differentiate between disk formation scenarios. Furthermore, 
the features that are present during the collapse are more eas¬ 
ily observed in the 450 pm continuum images than at longer 
wavelengths. It is difficult to test disk formation models to¬ 
ward highly inclined systems using continuum data since 
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both RSD and pseudo-disk formation show similar elongated 
features. 

The kinematical structures as revealed by the moment one 
maps of the synthetic molecular lines show a coherent blue- 
to red-shifted velocity gradient for both RSD models and for 
the pseudo-disk in the inner ~ 300 AU. However, the pseudo¬ 
disk shows a velocity gradient in north-south direction while 
the RSD shows an east-west gradient similar to the orien¬ 
tation of the flattened structure in the model. Moreover, the 
RSD formation in both 3D and 2D exhibits skewness in their 
moment one maps caused by the infalling rotating envelope 
component on larger scales which is absent in the case of 
pseudo-disk. The velocity gradient in the case of the pseudo¬ 
disk is a nearly straight line from the star to the large-scale 
structure since it is tracing the streams of material directly 
from the envelope onto the star. Thus, one can readily mis¬ 
take a pseudo-disk with an RSD unless one performs addi¬ 
tional analysis such as the peak position diagrams. 

Position-velocity (PV) diagrams constructed from and 
image cubes predict rotational signatures in both 
pseudo-disk and RSD formation, seen most prominently in 
data. This is due to the strength of rotation in the in¬ 
ner regions for both simulations. A combination of and 
lines is required to disentangle the RSD from a pseudo¬ 
disk. The signatures of infalling material in the pseudo-disk 
simulation are stronger in the lines. The velocity struc¬ 
ture constructed from the peak positions (peak-PV diagrams) 
are used to differentiate between the pseudo-disk and the 
RSD. Velocity structures described by i; cx and v oc r~^ 
are present in both cases of RSD formation whereas a flatter 
velocity profile is seen in the pseudo-disk case. We And that 
this conclusion is robust for different inclinations and rota¬ 
tions. 

The image cubes are convolved with large beams (> 9") to 
simulate single-dish observations probing the large-scale en¬ 
velope. The 2-1 and 3-2 line widths are similar be¬ 
tween the three simulations with FWHM ~ 1 km s“^or 
Doppler h of 0.4 km s“^ This is due to the fact that the 
emitting regions in the large-scale envelopes are similar. The 
observed FWHM values are larger than in those predicted 
from the simulations by a factor of 2. This suggests that the 
large-scale envelopes of low-mass embedded YSOs are sig¬ 
nificantly more turbulent than these models, which may be 
due to an interaction with a fast outflow (San Jose-Garcia 
et al. 2013). 

The comparison of the high-/ lines (^^CO 6-5, 10- 

9, and 9-8) indicates that the current simulations 

with RSD formation can reproduce the observed line widths 
solely due to the rotation -l- infall motions. On the other hand, 
the predicted line widths from the pseudo-disk are signifi¬ 
cantly smaller than the observed values. Thus, the mecha- 
nism(s) that are responsible for broadening the 6-5 and 9-8 
lines depend on the whether or not an RSD is present. If no 
RSD is present, the observations would imply an increas¬ 
ing level of turbulence with decreasing radii; if an RSD is 
present, turbulence would not be needed, although it can still 
be present in the disk at some level. 

Masses derived from continuum and molecular lines from 
simulations and analyzed in the same way as observations 
depend on the physical structure at small-scales. The disk 
masses obtained from the continuum flux in small beams (a 
few ") are generally in agreement with the true disk mass. 
Disk masses obtained from the integrated molecular line flux 
depend strongly on the physical structure. If the disk is small 


and cold with respect to the flattened inner envelope, inferred 
masses can be smaller by an order of magnitude from the 
true disk mass since the disk does not contribute to the in¬ 
tegrated line flux density. However, if the disk is large and 
warm enough to prevent significant freeze-out, the mass ob¬ 
tained from the low-/ lines can give a good estimate of the 
true RSD mass provided that the system is inclined toward 
us. Both the stellar masses and disk-to-envelope mass ratios 
are within a factor of 2 of the true masses. The presence of 
the RSD cannot be determined solely from continuum data 
alone, however the continuum flux provides a good estimate 
on the mass at small scales. Multiple spectrally and spatially 
resolved molecular line observations are needed to confirm 
the presence of RSD. 
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